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We  present  calculations  of  energetic,  electronic,  and  vibrational  properties  of  silicon  using  a  nonorthogonal 
tight-binding  (TB)  model  derived  to  fit  accurately  first-principles  calculations.  Although  it  was  fit  only  to  a  few 
high-symmetry  bulk  structures,  the  model  can  be  successfully  used  to  compute  the  energies  and  structures  of 
a  wide  range  of  configurations.  These  include  phonon  frequencies  at  high-symmetry  points,  bulk  point  defects 
such  as  vacancies  and  interstitials,  and  surface  reconstructions.  The  TB  parametrization  reproduces  experimen¬ 
tal  measurements  and  ab  initio  calculations  well,  indicating  that  it  describes  faithfully  the  underlying  physics 
of  bonding  in  silicon.  We  apply  this  model  to  the  study  of  finite  temperature  vibrational  properties  of  crystal¬ 
line  silicon  and  the  electronic  structure  of  amorphous  systems  that  are  too  large  to  be  practically  simulated  with 
ab  initio  methods. 


I.  INTRODUCTION 

As  the  capabilities  of  materials  simulations  increase,  so 
does  the  demand  for  methodologies  that  can  capture  the  im¬ 
portant  physics  accurately  while  being  fast  enough  to  simu¬ 
late  large  systems  for  long  periods  of  time.  Semiconductors 
such  as  silicon,  where  the  quantum-mechanical  nature  of  the 
electrons  mediating  the  directional  interatomic  bonding  is 
important,  have  been  especially  challenging  to  describe  ac¬ 
curately  using  empirical  potential  interactions.1  The  simplest 
method  that  captures  the  quantum-mechanical  nature  of  the 
electrons,  the  minimal  basis  tight-binding  (TB)  model,  is  be¬ 
coming  a  popular  method  for  simulating  such  systems. 2-10 
The  main  challenge  in  developing  these  models  has  been  the 
determination  of  the  Hamiltonian  matrix  elements  (and  in  the 
case  of  nonorthogonal  bases,  the  overlap  matrix  elements)  as 
a  function  of  interatomic  distance.  The  most  common  ap¬ 
proach  is  to  fit  the  results  of  either  total-energy  calculations 
or  band-structure  features  to  experiment  or  ab  initio  calcula¬ 
tions,  assuming  a  particular  functional  form  with  some  free 
parameters  for  the  distance  dependence  of  the  matrix  ele¬ 
ments.  Most  of  the  models  in  the  literature,  however,  suffer 
from  certain  shortcomings.  Many  models  assign  a  large  part 
of  the  total  energy  of  the  system  to  a  repulsive  pair  potential 
to  compensate  for  adopting  an  orthogonal  set  of  basis 
orbitals.2'5'6  Some  models  do  not  give  an  accurate  description 
of  the  band  stmcture  of  the  ground  state,7  or  the  energetics  of 
important  features  of  bulk  semiconducting  systems,  such  as 
point  defects.3,4  Some  models  were  intended  for  use  for  a 
specific  application,  and  were  optimized  and  tested  only  for 
geometries  relevant  to  that  application.9  Finally,  other  mod¬ 
els  use  a  very  large  number  of  free  parameters  and  fit  a  very 
large  number  of  structures,  leading  to  potentially  unphysical 
values  for  the  parameters  and  the  suspicion  that  the  good  fit 
to  the  data  set  may  not  guarantee  reliability.8 
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In  the  past,  TB  models  have  mainly  been  used  as  interpo- 
lative  schemes  assuming  that  the  configurations  that  were 
being  fit  will  encompass  the  configurational  space  where  the 
model  will  be  used.  In  this  paper  we  use  the  (nonorhogonal) 
NRL-TB  method,11  an  extrapolative  method  that  uses  param¬ 
eters  obtained  from  fitting  ab  initio  calculations  of  a  few 
high-symmetry  structures,  to  compute  the  energies  of  a  wide 
range  of  geometries  for  silicon.  We  find  that  the  results  com¬ 
pare  very  well  to  ab  initio  calculations  for  configurations  that 
are  substantially  different  from  those  included  in  the  fitting 
data  set.  This  increases  our  confidence  that  the  reason  for  the 
accurate  results  in  the  tested  configurations  is  that  the  phys¬ 
ics  underlying  the  model  is  correct. 

This  paper  is  organized  as  follows;  In  Sec.  II  we  describe 
the  functional  form  of  our  TB  parametrization  and  the  fitting 
data  set.  In  Sec.  Ill  we  discuss  applications  of  the  TB  model 
to  bulk  properties  such  as  the  diamond  lattice  band  structure 
and  energies  of  other  lattices,  point  defect  properties  such  as 
formation  and  relaxation  energies,  and  surface  energies  and 
reconstructions.  We  also  report  on  two  applications  of  the 
TB  model  that  would  be  impractical  with  ab  initio  methods, 
one  using  molecular  dynamics  (MD)  simulations  to  compute 
finite  temperature  properties  of  the  bulk  crystal,  and  another 
using  the  model  to  compute  electronic  properties  of  large 
amorphous  systems.  In  the  final  section  we  summarize  the 
results. 

II.  FUNCTIONAL  FORM  AND  FITTING 

In  this  paper  we  present  results  for  two  parametrizations, 
one  using  a  sp3  basis,  which  has  already  been  presented  in 
some  detail,10  and  another  using  a  sp^d5  basis.  Since  their 
functional  forms  are  nearly  identical  and  have  already  been 
presented,  we  will  only  give  a  brief  summary  here.  The  total 
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energy  of  the  system  is  written  as  the  sum  of  the  energies  of 
the  occupied  electronic  eigenstates.  The  onsite  Hamiltonian 
matrix  elements  vary  with  the  local  density,  allowing  the 
NRL-TB  method  be  fit  to  linearized  augmented  plane  wave 
(LAPW)  eigenvalues  that  have  been  shifted  so  that  the 
LAPW  total  energy  is  equal  to  the  eigenvalue  sum.  There¬ 
fore  all  of  the  contributions  to  the  total  energy  are  accounted 
for  in  the  eigenvalue  sum  and  the  addition  of  a  repulsive  pair 
potential,  a  feature  common  to  most  TB  models,  is  not 
needed. 

The  energies  of  the  electronic  states  and  the  correspond¬ 
ing  eigenvectors  are  the  solutions  to  a  generalized  eigenvalue 
equation  with  Hamiltonian  and  overlap  matrix  elements  pa¬ 
rametrized  as  follows:  The  basis  used  to  describe  the  Hamil¬ 
tonian  and  overlap  matrices  is  a  set  of  one  s,  three  p,  and  for 
one  of  the  parametrizations  five  d  orbitals  around  each  atom, 
with  all  interactions  assumed  to  be  in  the  two-center 
approximation.12  A  local  atomic  density  at  atom  i  is  defined 
as 

P,  =  2  e-^r^fURj-R'l),  (1) 

j 


TABLE  I.  Parameters  for  the  sp3  basis  tight-binding  model. 


X 

Orbital 

Onsite  parameters 

1.1036 

a  (Ry)  ft  (Ry)  y  (Ry) 

x  (Ry) 

s 

-0.0532 

-0.9076 

-8.8308 

56.5661 

p 

0.3579 

0.3036 

7.0922 

-77.4786 

Hamiltonian  matrix  parameters 

Interaction 

a  (Ry) 

b  (Ry/a.u.) 

c  (Ry/a.u.2) 

g  (a.uf1/2) 

I_l 

11  sscr 

219.5608 

-16.2132 

-15.5049 

1.2644 

i_i 

11  spa 

10.1279 

-4.4039 

0.2267 

0.9227 

i_i 

11  ppa 

-22.9590 

1.7208 

1.4191 

1.0314 

i_i 

‘  ‘  pp  TT 

10.2654 

4.6718 

-2.2162 

1.1113 

Overlap  matrix  parameters 

Interaction 

t  (a.u.  1 ) 

q  (a.u.  2) 

r  (a.u.  3) 

u  (a.u71/2) 

c 

u  ssa 

5.1576 

0.6600 

-0.0815 

1.1081 

C 

u  spa 

8.8736 

-16.2408 

5.1823 

1.2407 

C 

°  PPa 

11.2505 

-1.1701 

-1.0591 

1.1376 

C 

u  PPTT 

-692.1842 

396.1532 

-13.8172 

1.5725 

where  R ,  is  the  position  of  atom  i  and  X  is  a  fitting  param¬ 
eter.  The  cutoff  function  f(R)  is  given  by 


f(R)  = 


1  +exp 


R  —  R ..  +  5L. 


-l 


R^R, 


(2) 


0  R>RC, 


where  Rc  is  12.5  a.u.  and  Lc  is  0.5  a.u.  The  onsite  matrix 
elements  are  given  in  terms  of  the  local  atomic  density  p,  as 


hn~ai  +  Pipf3+7ipf3+XiPh  (3) 

where  /  is  the  orbital  type  index  (s,  p,  or  d)  and  at ,  f3t ,  y; , 
and  Xi  are  fitting  parameters.  The  distance  dependent  parts  of 
the  two-center  Hamiltonian  matrix  elements  are  given  by 

Hii'M.(R)  =  (ai,t/l+bll,fJ,R  +  Cn'/J,R2)exp(-gl,/iR)f(R), 

(4) 

where  /  and  / '  are  orbital  type  indices,  /j.  is  an  index  for  the 
type  of  interaction  between  orbitals  (cr,  it,  or  S),  and  the 
parameters  aiv ^  bp,^,  Cjp^,  and  are  fitting  param¬ 
eters.  The  distance  dependent  parts  of  the  overlap  matrix 
elements  are 


Sw  n(R)-{8ir+tii'fJ.  +  qil'  f±R  +  ru,  ^R1) 

Xexp(-M2I(tR )/(/?),  (5) 

where  tpr^,  qw^,  rwp,  and  uni ^  are  fitting  parameters, 
and  Sip  is  the  Kronecker  delta.  Note  that  the  symbols  for 
some  of  the  parameters  are  different  from  those  used  in  Ref. 
10.  The  overlap  matrix  elements  have  similar  functional 
form  to  the  Hamiltonian  matrix  elements,  but  are  constrained 
to  go  to  the  correct  value,  zero  or  one,  at  zero  interatomic 
separation.  The  angular  dependence  of  the  Hamiltonian  and 
overlap  matrix  elements  is  the  standard  two-center  Slater- 
Koster  form.12 

The  41  parameters  used  by  the  functional  form  for  the  sp 3 
basis  parametrization  are  fit  to  four  high-symmetry  crystal 


lattices:  simple  cubic  (sc),  face-centered  cubic  (fee),  body- 
centered  cubic  (bcc),  and  the  diamond  structure.  The  fitting 
data  set  includes  both  the  total  energy  and  band  structure  of 
each  lattice,  as  computed  by  LAPW  ab  initio  density- 
functional  theory  (DFT)  calculations  in  the  local-density  ap¬ 
proximation  (LDA)  for  a  wide  range  of  volumes  around  the 
energy  minimum.  The  diamond  lattice  data  included  the  wid¬ 
est  range  of  volumes,  from  12.2  A3/atom  to 
22.7  A3/atom.  The  sc  lattice  structures  ranged  from 

12.6  A3/atom  to  18.5  A3/atom,  the  fee  lattice  from 

12.7  A  3/ atom  to  15.0  A  3/ atom,  and  the  bcc  lattice  from 
13.0  A3/atom  to  16.0  A3/atom.  The  best-fit  root-mean- 
square  (rms)  error  of  the  valence-band  energies  for  the  dia¬ 
mond  lattice  is  0.12  eV,  and  the  rms  error  for  the  crystal 
lattice  total  energies  is  0.020  eV. 

The  sp3d5  basis  parametrization  has  69  parameters, 
which  were  fit  to  the  diamond  lattice  band  structure  at  vol¬ 
umes  ranging  from  13.5  A3/atom  to  22.7  A3/atom.  This 
set  of  parameters  does  not  allow  for  any  Hamiltonian  or 
overlap  matrix  elements  between  different  d  orbitals,  but  al¬ 
lows  all  interactions  between  s  and  p  orbitals,  as  well  as  s 
—  d  and p—d  interactions.  Since  this  model  is  optimized  for 
accuracy  in  the  band  structure,  we  adjusted  the  DFT/LDA 
calculations,  which  predict  an  indirect  gap  of  0.5  eV,  by 
applying  a  uniform  shift  of  0.67  eV  to  the  conduction  bands 
of  the  ideal  volume  diamond  lattice,  matching  the  band  gap 
to  the  experimental  result.1314  For  the  other  lattice  constants 
in  the  fit  we  shifted  the  conduction  band  so  that  the  gap  was 
scaled  up  by  a  factor  of  1.17/0.5  =  2.34.  The  shift  amount 
increased  monotonically  from  0.67  eV  for  larger  volumes 
and  decreased  monotonically  for  smaller  volumes.  The  best 
fit  rms  error  for  the  diamond  lattice  valence  band  and  lower 
two  conduction  bands  is  0.21  eV,  and  the  rms  error  for  the 
diamond  lattice  total  energies  is  0.004  eV. 

The  parameters  that  result  from  this  fit  for  the  sp3  and 
sp3d 5  basis  models  are  listed  in  Tables  I  and  II,  respectively. 
The  sp 3  basis  Hamiltonian  and  overlap  matrix  elements  are 
plotted  in  Fig.  1,  and  the  onsite  matrix  elements  for  the  dia- 
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TABLE  II.  Parameters  for  the  sp3d 5  basis  tight-binding  model. 


Onsite  parameters 

\ 

1.1108 

Orbital 

a  (Ry) 

yS  (Ry) 

r  (Ry) 

x  (Ry) 

s 

-0.0555 

-1.1131 

-7.3201 

74.8905 

p 

0.4127 

-0.0907 

5.3155 

-44.0417 

d 

0.9691 

-0.9151 

-5.9743 

602.0289 

Hamiltonian  matrix  parameters 

Interaction 

a  (Ry) 

b  (Ry/a.u.)  c 

(Ry/a.u.2) 

g(a.u.-'/2) 

LJ 

11  sscr 

234.6937 

-18.6013 

-15.0266 

1.2502 

LJ 

11  spa 

9.5555 

-4.1279 

0.2499 

0.8761 

LJ 

n  ppa 

-22.6782 

1.3611 

1.3879 

1.01655 

H pprr 

-1.5942 

4.7914 

-1.5693 

1.1030 

H  sd  a 

-7571.4416 

2.2354 

7.0122 

1.6234 

H  pda 

-1.8087 

-3.4695 

-7.7637 

1.6294 

H pdrr 

0.8933 

0.1058 

-0.0224 

0.8217 

Overlap  matrix  parameters 

Interaction 

t  (a.uf1) 

q  (a.u.  2)  r 

■  (a.u.-3) 

xu  (a.u.  I/2) 

2.4394 

0.9091 

-0.0749 

1.0590 

C 

^ spa 

-12.0027 

-14.6860 

6.1856 

1.2218 

C 

^ ppa 

13.9608 

-1.1961 

-1.2606 

1.1118 

C 

'-’pprr 

188.0012 

-143.3625 

33.5043 

1.4340 

$  sda 

11.4724 

-0.4454 

-0.5838 

1.0598 

S  pda 

-0.6071 

0.05789 

0.0221 

0.8130 

^  pdrr 

-2.1340 

-0.5209 

-0.0948 

1.0580 

mond  structure  are  plotted  in  Fig.  2.  The  variation  of  the 
onsite  matrix  elements  with  nearest-neighbor  distance  is 
structure  dependent  because  they  have  a  nonlinear  depen¬ 
dence  on  the  density,  which  is  itself  a  structure-dependent 
quantity.  Note  that  the  sscr,  ppa,  and  pp  it  Hamiltonian  and 
overlap  parameters  have  the  expected  sign,  while  the  spa 
parameters  are  opposite  in  sign  to  the  usual  convention.12 
However,  this  sign  is  not  physically  meaningful,  since  it  is 
determined  by  the  (arbitrary)  choice  of  sign  of  the  s  and  p 
basis  orbitals,  and  does  not  affect  the  eigenvalues  or  energies 
computed  with  the  model.  Cohen  et  al.  have  also  used  a 
similar  method  to  generate  parameters  for  silicon.9  They 
used  a  somewhat  different  functional  form  and  concentrated 
their  fitting  and  tests  on  high-pressure  phases.  Since  we  are 
interested  in  applying  this  TB  model  to  complex  structures, 
including  defects  and  surfaces,  but  at  or  near  atmospheric 
pressure,  we  have  employed  a  different  set  of  geometries  for 
our  fitting  data  set. 


III.  APPLICATIONS 

To  test  the  transferability  of  the  sp3  parameters  we  com¬ 
puted  the  total  energy  of  a  range  of  structures  important  for 
condensed  phases  of  silicon  including  bulk  systems,  point 
defects,  and  surfaces.  First  we  review  the  diamond  lattice 
band  structure,  cohesive  energies  of  a  number  of  bulk  lattices 
as  a  function  of  volume,  and  the  elastic  constants  of  the 
diamond  structure,  as  presented  in  Ref.  10.  To  address  the 
issue  of  improving  the  diamond  lattice  band  structure  we 


r  (Angstrom) 


r  (Angstrom) 


FIG.  1.  Hamiltonian  matrix  elements  (upper  panel)  and  overlap 
matrix  elements  (lower  panel)  for  the  sp3  parametrization  plotted  as 
a  function  of  interatomic  distance. 


present  an  electronic  structure  calculation  using  the  sp3d5 
parametrization.  We  expand  our  analysis  of  the  sp3  param¬ 
eters  to  include  phonon  spectra  at  several  high-symmetry 
points.  As  a  more  stringent  test  we  use  this  parametrization 
to  compute  the  energies  of  some  lower-symmetry  configura¬ 
tions.  For  the  bulk  we  simulate  important  point  defects,  in¬ 
cluding  several  low-energy  interstitial  configurations,  the  va¬ 
cancy,  and  the  concerted  exchange  pathway  for  diffusion. 
For  the  (100)  and  (111)  surfaces  we  compute  ideal  surface 
energies  as  well  as  relaxation  energies  for  a  number  of  re¬ 
constructions.  From  MD  simulations  we  compute  the  mean- 
squared  atomic  displacement  for  a  range  of  temperatures,  the 
vibrational  density-of-states,  and  the  phonon-dispersion  rela¬ 
tions.  Finally,  we  use  the  sp3d 5  basis  model  to  study  the 
electronic  structure  of  large  amorphous  systems. 

A.  Bulk 

The  band  structure  of  the  diamond  structure  lattice,  which 
was  part  of  the  fitting  data  set,  is  shown  in  Fig.  3.  The  va¬ 
lence  band  is  in  very  good  agreement  with  LAPW  calcula¬ 
tions.  The  conduction  band  is  not  as  well  described,  with  the 


(Angstrom) 


FIG.  2.  Onsite  matrix  elements  for  the  sp 3  parametrization  for 
the  diamond  structure  plotted  as  a  function  of  the  lattice  constant. 
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r  X  W  L  r  K 


energy  (eV) 


FIG.  3.  Band  structure  of  Si  in  the  diamond  lattice  computed 
with  the  sp 3  (upper  panel)  and  sp3d5  (middle  panel)  bases,  and  the 
density-of-states  for  the  sp3d5  model  (bottom  panel).  Dashed  lines 
are  TB  results,  solid  lines  are  DFT/LDA  results  (a  rigid  shift  has 
been  applied  to  the  DFT/LDA  conduction-band  results  used  for  the 
sp3d5  basis  model).  All  energies  are  referred  to  as  the  valence-band 
maximum. 

minimum  indirect  gap  of  1.02  eV  appearing  at  the  L  point 
rather  than  at  approximately  three  fourths  of  the  way  from 
the  r  to  the  X  point  as  first-principles  calculations  and  ex¬ 
periments  have  established.  The  size  of  the  gap  is  somewhat 
smaller  than  the  experimental  value  of  1.17  eV,15  although  it 
is  larger  than  the  DFT/LDA  prediction  to  which  it  was  fitted. 

We  have  addressed  the  issue  of  obtaining  a  better  fit  of 
the  gap  and  the  conduction  band  and  came  to  the  following 
conclusions.  The  addition  of  five  d  orbitals  to  the  basis  im¬ 
proves  the  fit  of  the  conduction  band,  as  can  be  seen  in  Fig. 
3.  The  lowest-energy  conduction  band  is  nearly  perfect  when 
compared  to  a  LAPW  calculation  with  a  rigid  0.67  eV  shift 
of  the  conduction  bands,1314  and  the  higher  bands  are  also 
closer  to  the  LAPW  calculation  than  with  the  sp 3  basis 
model.  The  valence  bands  are  very  well  described,  although 
the  lowest  band  at  the  I  point  is  too  fiat.  The  density-of- 
states  (DOS),  including  its  decomposition  into  contributions 
from  different  angular  momentum  states  (which  is  found  us¬ 
ing  TB  eigenvectors  that  were  not  fitted),  is  also  in  very  good 
agreement  with  DFT/LDA  calculations.  The  three  peaks  in 
the  valence  band  are  clear,  as  is  the  decomposition  into 


FIG.  4.  Total  energy  vs  volume  for  the  diamond  structure  as 
well  as  a  number  of  crystal  structures  for  Si  not  included  in  the 
fitting,  computed  using  the  TB  model.  Si34  is  a  clathrate  structure, 
BC8  is  the  body  centered  eight-atom  structure,  yS-Sn  is  the  structure 
of  the  /3  phase  of  tin,  and  hep  is  the  hexagonal-close  packed  struc¬ 
ture. 

mainly  s  character  in  the  lowest  peak,  mixed  s  and  p  charac¬ 
ter  in  the  middle  peak,  and  p  character  in  the  third  peak. 
There  is  very  little  d  character  in  the  valence  band,  while  the 
conduction  band  is  mainly  of  mixed  p  and  d  character,  with 
smaller  ,v  contribution. 

To  obtain  such  a  good  fit  for  the  band  structure  the  sp3d5 
model  was  fit  only  to  the  full  diamond  lattice  band  structure 
at  all  volumes.  The  lack  of  other  structures  and  energy  infor¬ 
mation  in  the  fitting  data  set  deteriorates  the  energetics  of  the 
model.  We  decided  that  the  best  compromise  is  to  use  the 
minimal  sp3  basis  for  all  total-energy  calculations  presented 
in  this  paper,  and  to  use  the  sp3d5  to  compute  the  electronic 
structure  of  amorphous  silicon  presented  later  in  this  section. 
We  note  in  passing  that  in  his  book,  Papaconstantopoulos16 
was  able  to  obtain  a  good  fit  of  the  conduction  band  near  the 
gap  with  a  sp 3  basis  model.  However,  that  work  differs  from 
the  present  approach  in  two  important  ways:  first,  it  utilizes 
three-center  integrals  and  second,  it  treats  the  Hamiltonian 
and  overlap  matrix  elements  for  the  first  three  neighbor 
shells  as  independent  parameters,  rather  than  giving  them  as 
an  analytical  functional  form  that  varies  with  distance.  These 
differences  provide  the  flexibility  that  produces  a  better  fit  to 
the  conduction  states. 

The  total  energies  as  a  function  of  volume  for  a  wide 
range  of  structures  are  shown  in  Fig.  4,  and  their  equilibrium 
structural  and  energy  properties  are  listed  in  Table  III.  All  of 
the  structures  have  higher  energy  than  the  diamond  structure, 
including  some  low-energy,  rarely  examined  theoretical 
phases  such  as  hexagonal  diamond  and  the  clathrate 
structures.17  The  TB  model  reproduces  the  LAPW  results 
very  well  for  the  four  structures  to  which  it  was  fit,  as  can  be 
seen  from  the  equilibrium  energies,  volumes,  and  bulk 
moduli  listed  in  Table  IV.  These  quantities  were  computed 
using  a  Birch  fit18  to  the  fitting  data  and  the  sp3  TB  model 
calculations.  The  elastic  constants  of  the  ground-state  dia¬ 
mond  structure  are  listed  in  Table  V.  Those  that  do  not  in¬ 
volve  shear,  cu  and  c12,  are  within  22%  of  LAPW  calcula¬ 
tions  and  14%  of  experiment.19  The  shear  elastic  constant 
computed  without  allowing  for  relaxation  of  the  internal  de¬ 
grees  of  freedom  of  the  two-atom  unit  cell  c*4,  is  34%  larger 
than  the  LAPW  result.  Allowing  for  the  internal  relaxation 
brings  c44  within  19%  of  experiment. 19  A  detailed  compari- 
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TABLE  III.  Equilibrium  energies  and  structural  features  of  hy¬ 
pothetical  crystal  lattices  for  Si  computed  with  the  sp 3  TB  model.  E 
is  equilibrium  energy  relative  to  the  diamond  structure  in  eV  per 
atom,  V  is  the  volume  in  A  3  per  atom,  and  da  is  the  unit  cell 
aspect  ratio.The  internal  structural  parameter  x  for  the  H-Dia  struc¬ 
ture  is  the  position  (y,f,x)  of  the  atom  at  site  4/,  and  for  the  BC8 
structure  is  the  position  (x,x,x)  of  the  atom  at  site  16c.  Notation  for 
the  lattice  types  is  the  same  as  in  Fig.  4 


Lattice 

E 

V 

c/a 

X 

Dia 

0.000 

19.97 

H-Dia 

0.021 

19.94 

1.647 

0.0630 

BC8 

0.229 

18.05 

0.1008 

/3-Sn 

0.357 

14.56 

0.5278 

SH 

0.389 

14.93 

0.9479 

CO 

0.480 

13.77 

0.5917 

HCP 

0.498 

14.04 

1.637 

son  of  the  energetic  and  structural  parameters  of  two  clath- 
rate  structures,  Si34  and  Si4f) ,  comparing  results  of  our  TB 
model  with  experiment,  DFT  calculations,  and  results  of  an 
orthogonal  TB  model,  is  shown  in  Table  VI.  The  energy 
differences  between  the  clathrates  and  the  diamond  structure 
are  lower  than  in  DFT  calculations,  although  they  are  of  the 
correct  sign.  The  structures,  including  both  the  lattice  con¬ 
stant  and  the  internal  structural  parameters  of  the  basis,  are 
within  1%  of  the  experimental  values.20'21  This  agreement  is 
as  good  as  that  provided  by  DFT/LDA  plane-wave 
calculations22  or  by  an  ab  initio  localized  orbital  method,17 
and  substantially  better  than  the  orthogonal  TB  model  of 
Goodwin  et  al.2  tested  by  Kahn  and  Lu.23 

Phonon  frequencies  at  high-symmetry  points  in  the  Bril- 
louin  zone  (BZ)  computed  with  the  TB  model  using  the  fro¬ 
zen  phonon  approximation  are  compared  with  experimen¬ 
tally  measured  values  in  Table  VII.24  The  agreement  is  quite 
good:  the  TB  results  are  within  15%  of  experiment  for  all  but 
three  of  the  modes,  the  X3,  L , ,  and  W2  modes,  which  are  off 
by  about  25%,  30%,  and  60%,  respectively.  While  this  good 
description  of  the  phonon  spectra  is  a  nontrivial  test  of  the 

TABLE  IV.  Equilibrium  energy  relative  to  the  diamond  struc¬ 
ture  ( E )  in  eV  per  atom,  volume  ( V )  in  A  3  per  atom,  and  bulk 
modulus  (B)  in  GPa,  for  the  lattice  structures  in  the  fitting  data  set, 
computed  with  the  sp2  TB  model  and  with  LAPW  DFT/LDA. 


Lattice 

TB 

DFT/LDA 

Diamond 

E 

0.000 

0.000 

V 

19.97 

19.78 

B 

108.3 

96.4 

SC 

E 

0.279 

0.338 

V 

15.17 

15.76 

B 

101.5 

105.6 

fee 

E 

0.495 

0.484 

V 

14.28 

13.85 

B 

117.1 

93.54 

bcc 

E 

0.474 

0.439 

V 

13.58 

14.08 

B 

88.56 

111.3 

TABLE  V.  Elastic  constants  in  GPa  for  the  diamond  structure, 
computed  with  the  sp3  TB  model,  LAPW  DFT/LDA  calculations, 
and  experiment.  c44  is  the  shear  elastic  modulus  computed  without 
allowing  for  relaxation  of  the  internal  degrees  of  freedom  in  the 
two-atom  unit  cell. 


TB 

DFT/LDA 

Expt. a 

C11 

179 

152 

166 

c  12 

73 

60 

64 

r* 

c44 

135 

101 

c44 

95 

80 

“From  Ref.  19. 


model,  it  represents  only  infinitesimal  deviations  of  atomic 
positions  from  the  diamond  structure,  which  was  part  of  the 
fitting  data  set.  In  the  following  sections  we  show  that  this 
TB  model  can  also  accurately  describe  the  energies  of  con¬ 
figurations  that  are  substantially  different  from  those  in  the 
fitting  data  set. 

B.  Point  defects 

The  ground-state  structure  for  silicon  is  the  covalently 
bonded,  open  network  of  the  diamond  structure.  Since  strong 
covalent  bonds  in  the  ideal  lattice  allow  for  little  atomic 
motion  at  temperatures  below  the  melting  point,  processes 
such  as  diffusion  are  dominated  by  point  defects,  which  are 
far  more  mobile  than  perfectly  bonded  atoms.25  The  forma¬ 
tion  energy  of  such  defects  strongly  influences  their  concen¬ 
trations  and  is  therefore  an  important  material  property.  In 

TABLE  VI.  Energy  differences  A E  relative  to  the  diamond 
structure  and  structural  parameters  for  the  two  optimized  clathrate 
structures,  Si34  and  Si46 .  A E  is  given  in  eV/atom,  the  lattice  con¬ 
stant  a0  is  given  in  A  ,  and  the  internal  parameters  are  given  in 
terms  of  the  lattice  constant.  Using  the  notation  of  Ref.  17,  the 
parameters  of  the  Si34  structure  are  the  position  ( xe  ,xe  ,xe)  of  the 
atom  at  site  32a,  and  the  position  (xg ,xg  ,zg)  of  the  atom  at  site 
96 g.  The  parameters  of  the  Si46  structure  are  the  position  (xt  ,xt  ,x,) 
of  the  atom  at  site  16/  and  the  position  (0,yk  ,zk )  of  the  atom  at  site 
2Ak.  PW  denotes  plane-wave  basis  DFT/LDA  calculations  from 
Ref.  22,  LO  denotes  local  orbital  DFT/LDA  calculations  from  Ref. 
17,  OTB  denotes  the  orthogonal  tight-binding  results  from  Ref.  23, 
and  NRL-TB  denotes  the  present  paper.  Note  that  the  experimental 
samples  are  actually  of  NavSi34  with  x<  11,  and  Na8Si46 . 


Exper. 

PW 

LO 

OTB 

NRL-TB 

Si34 

A  E 

0.077 

0.055 

0.004 

14.62 

14.55 

14.864 

14.479 

14.543 

0.219 

0.2171 

0.2174 

0.2020 

0.2171 

xt 

0.183 

0.1825 

0.1824 

0.1823 

0.1824 

0.371 

0.3705 

0.3701 

0.3703 

0.3704 

St48 

A  E 

0.069 

0.018 

ao 

10.19 

10.355 

10.055 

10.089 

xt 

0.183 

0.1837 

0.1835 

0.1835 

yk 

0.116 

0.1172 

0.1174 

0.1171 

Zk 

0.310 

0.3077 

0.3077 

0.3082 
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TABLE  VII.  Phonon  frequencies  (in  cm-1)  at  high-symmetry 
points  of  the  BZ  computed  with  the  TB  model  and  measured  ex¬ 
perimentally.  A  =  (0,0,77 /a)  and  2  =  (0,ir/2a,T7/2a). 


TB 

Exp.  a 

r 

531 

518 

Xi  (L) 

405 

415 

V3  (L) 

508 

406 

v4  (T) 

160 

150 

Li  (L) 

553 

417 

L2  (L) 

333 

383 

t3+  (T) 

533 

487 

l3-  (T) 

127 

114 

IV, 

371 

403 

W2 

221 

140 

wv 

514 

470 

A,  (L) 

235 

237 

^3 

518 

500 

a4 

144 

131 

a4, 

525 

474 

238 

270 

2,, 

528 

500 

s2 

150 

141 

s3 

210 

201 

Sy 

452 

464 

24 

525 

480 

aFrom  Ref.  24. 


Table  VIII  we  list  the  formation  energies  of  three  interstitial 
configurations,  the  tetrahedral,  hexagonal,  and  (110)  split, 
and  of  the  vacancy.  All  defect  energies  were  computed  using 
a  16.29  A  cube  cell  with  216±  1  atoms,  and  sampling  the  BZ 
at  the  T  point.  To  compute  the  relaxed  configurations  a 
conjugate-gradient  algorithm  was  used,26  with  the  atomic  po¬ 
sitions  relaxed  until  the  force  on  each  atom  was  less  than  3 
meV/A  .  In  agreement  with  ah  initio  calculations  the  (110) 
split  is  the  lowest-energy  interstitial  configuration,  followed 
by  the  tetrahedral  and  hexagonal  configurations.27'28  The  for¬ 
mation  and  relaxation  energies  of  all  three  interstitial  con¬ 
figurations  are  within  10%  of  the  range  of  DFT/LDA 
calculations.27  29  The  formation  energy  of  the  ideal  vacancy 
is  also  accurately  predicted,  although  its  relaxation  energy  is 
about  twice  as  large  as  DFT/LDA  calculations  predict.27’30 

The  relaxed  geometries  are  in  approximate  agreement 
with  ah  initio  results,  ’  ’  although  there  are  some  differ- 

TABLE  VIII.  Formation  energies  E'feal  and  relaxation  energies 
A  Erf,ax  ,  in  eV,  for  point  defects  computed  using  the  TB  model  and 
comparison  to  DFT/LDA  results.  Since  the  structure  of  the  ideal 
split  interstitial  is  not  uniquely  defined,  the  energy  listed  under 
E'feat  is  actually  the  relaxed  formation  energy. 


TB  LDA a 


rrideal 

Ef 

A  Er/lax 

rrideal 

Ef 

A  Er/lax 

(110)  split  interstitial 

3.7 

3.3 

Tetrahedral  interstitial 

4.8 

0.3 

3.7-4. 8 

0.1 -0.2 

Hexagonal  interstitial 

5.4 

0.5 

4.3-5.0 

0.6-1. 1 

Vacancy 

4.2 

1.0 

3. 3-4. 3 

0.4-0.6 

“From  Refs.  27-30. 


FIG.  5.  Structure  of  the  tetrahedral  interstitial  in  ideal  (white 
spheres)  and  relaxed  (black  spheres)  configurations,  with  bonds 
connecting  atoms  within  2.5  A  .  The  geometry  is  viewed  along  the 
(110)  direction  with  the  (111)  direction  pointing  up. 

ences.  The  tetrahedral  interstitial  reduces  its  symmetry,  with 
three  of  the  four  atoms  that  surround  it  relaxing  outward  and 
the  fourth  relaxing  inward,  while  the  interstitial  atom  itself 
moves  parallel  to  the  fourth  atom.  A  diagram  of  this  configu¬ 
ration,  shown  in  Fig.  5,  makes  clear  that  the  relaxed  intersti¬ 
tial  atom  has  moved  part  way  (0.34  A)  towards  a  hexagonal 
site  without  significantly  stretching  (<2.5%)  any  of  its 
bonds.  This  geometry  differs  from  the  DFT/LDA  work 
where  the  relaxed  interstitial  had  the  full  symmetry  of  the 
initial  configuration,27  although  whether  those  calculations 
were  restricted  to  maintain  tetrahedral  symmetry  or  simply 
found  no  energy  gain  from  breaking  the  symmetry  is  not 
stated.  The  hexagonal  interstitial  retains  its  ideal  hexagonal 
symmetry  with  the  six  ring  atoms  moving  outward  by  0. 1  A, 
in  perfect  agreement  with  the  DFT/LDA  results.27  The  four 
atoms  around  vacancy  relax  inward  by  about  0.28  A  ,  in 
good  agreement  with  DFT/LDA  calculations,  which  find  an 
inward  relaxation  of  0.25  A  .  However  the  structure  keeps 
the  full  tetrahedral  symmetry  rather  than  reducing  to  the 
symmetry  of  the  tetragonal  structure  that  ah  initio  calcula¬ 
tions  predict.27,30 

To  test  the  accuracy  of  the  TB  model  in  describing  the 
breaking  of  bonds  within  a  relatively  undisturbed  solid  we 
computed  the  energy  for  the  concerted  exchange  pathway  for 
diffusion.3132  The  energy  along  the  path  for  the  ideal  and 
relaxed  configurations,  as  well  as  DFT/LDA  results,31,32  are 
plotted  in  Fig.  6.  The  ideal  energy  is  within  22%  of  the 
DFT/LDA  calculation  throughout  the  path,  with  no  spurious 
minima.  The  relaxation  energy  is  substantially  too  high,  al¬ 
though  the  relaxed  length  of  the  bond  between  the  diffusing 
atoms  of  2.15  A  is  nearly  identical  to  the  DFT/LDA  result. 

Point  defect  configurations  include  substantial  deviations 
from  the  ideal  lattice  geometry  and  several  inequivalent 
atomic  sites.  In  such  a  situation  it  is  possible  for  charge  to  be 
transfered  between  atoms.  If  this  charge  transfer  is  substan¬ 
tial,  the  applicability  of  a  model  with  no  Coulombic  interac¬ 
tion  or  charge  self-consistency  may  be  in  doubt.  The  varia¬ 
tion  of  the  onsite  energies  in  our  TB  model  could  potentially 
exacerbate  this  effect.  To  measure  the  amount  of  charge 
transfer  we  performed  a  Mulliken  population  analysis  on  the 
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Angle  (deg) 

FIG.  6.  Energy  along  the  concerted  exchange  pathway  for  the 
diffusion  of  atoms  without  vacancies  or  interstitials.  DFT/LDA  cal¬ 
culations  of  the  ideal  configuration  from  Ref.  31  are  plotted  in  a 
solid  line,  TB  calculations  of  the  ideal  configurations  in  a  dashed 
line,  and  TB  calculations  of  the  relaxed  configurations  in  a  dot-dash 
line.  The  relaxed  DFT/LDA  value  is  only  available  for  the  saddle 
point. 

vacancy,  the  tetrahedral  interstitial,  and  the  hexagonal  inter¬ 
stitial.  The  deviation  of  the  charge  from  neutrality  (four  elec¬ 
trons  per  atom)  was  modest,  less  than  half  an  electron  in 
every  case,  and  as  low  or  lower  for  the  relaxed  configura¬ 
tions  as  for  the  ideal  ones.  By  comparison,  a  nonorthogonal 
TB  model7  with  constant  onsite  energies  predicted  somewhat 
smaller  charge  transfers  (by  30%  to  50%),  except  for  the 
ideal  tetrahedral  interstitial  where  it  predicted  —  0.65<?  as 
compared  with  —  0.35e  on  the  interstitial  atom.  This  com¬ 
parison  indicates  that  the  variation  of  the  onsite  elements 
does  not  substantially  increase  the  charge  transfer,  which  are 
neglected  by  most  TB  models. 

C.  Surfaces 

The  properties  of  the  silicon  surface  are  critical  for  pro¬ 
cesses  such  as  surface  growth  and  etching.  Since  atoms  on  a 
surface  have  an  asymmetric  environment  and  lower  coordi¬ 
nation  than  in  the  bulk,  the  resulting  configurations  are  quali¬ 
tatively  different  from  the  types  of  geometries  that  this  TB 
parametrization  was  fitted  to.  Therefore,  surface  structures 
provide  an  important  test  of  the  transferability  of  the  param¬ 


eters.  We  compute  the  surface  energies  of  the  ideal  (100)  and 
(111)  surfaces,  as  well  as  the  relaxation  energies  for  some 
simple  but  representative  reconstructions  of  these  surfaces. 
For  both  surfaces  we  use  a  symmetric  slab  configuration  [24 
layers  for  the  (100)  surface,  6  bilayers  for  the  (111)  surface] 
and  a  set  of  k  points  equivalent  to  a  4  X  4  mesh  in  the  full 
planar  BZ  of  the  1X1  surface  unit  cell.  For  the  (100)  surface 
we  examined  the  2X1  buckled  dimer  reconstruction,  and  for 
the  (111)  surface  we  examined  2X2  reconstructions  with 
adatoms  at  three  inequivalent  sites,  the  7  4 ,  H3,  and  B2.  The 
results  of  these  calculations  are  listed  in  Table  IX. 

The  energetics  of  the  ideal  surfaces  are  in  agreement  with 
DFT/LDA  calculations,33’34  with  the  (111)  surface  about  0.8 
eV  lower  in  energy  than  the  (100)  surface.  Both  the  relax¬ 
ation  energy  and  the  structure  of  the  buckled  dimer  recon¬ 
struction  of  the  (100)  surface  are  in  good  agreement  with 
DFT/LDA  calculations.33’35’36  One  of  the  three  adatom  re¬ 
constructions  of  the  (111)  surface,  with  the  adatom  at  the  T4 
site,  is  related  to  the  dominant  features  of  the  ground-state 
7X7  reconstruction.37  The  //3  adatom  site  is  a  metastable 
position  that  is  involved  in  the  diffusion  of  adatoms.  The 
energy  of  the  B2  site,  which  lies  halfway  along  the  path 
connecting  the  T4  and  H2  sites,  determines  the  barrier  for 
migration  between  them. 

The  TB  model  predicts  the  same  ordering  in  energy  as 
DFT/LDA  calculations  for  the  three  adatom  sites,34’38  39  an 
improvement  over  previous  nonorthogonal  TB  models.7  The 
agreement  is  not  quantitative,  however,  as  the  TB  model 
overestimates  the  binding  energy  of  the  adatom  at  the  T 4  site 
by  about  50%.  In  this  position  the  adatom  forms  a  bond  of 
length  of  2.34  A  to  the  atom  underneath,  and  bonds  of  length 
2.42  A  to  the  three  next-nearest-neighbors.  DFT/LDA  calcu¬ 
lations  predict  corresponding  bond  lengths  of  2.43  A  and 
2.47  A  ,34  The  surface  atom  that  does  not  form  bonds  to  the 
adatom,  called  the  rest  atom,  makes  three  bonds  of  length 
2.38  A  with  angles  of  100.5°,  in  very  good  agreement  with 
the  DFT/LDA  result  of  2.34  A  long  bonds  with  angles  of 
99.9°.  The  binding  energy  of  the  adatom  in  the  H2  site  is 
somewhat  underestimated  as  compared  with  DFT/LDA 
calculations.34-39  In  this  site  the  adatom  makes  three  2.39  A 


TABLE  IX.  Surface  energies  y  (in  eV  per  1X1  surface  unit  cell),  for  the  ( 100)  and  (111)  surfaces  of  Si, 
relaxation  energies  A  y  (in  eV  per  1X1  surface  unit  cell)  relative  to  the  ideal  surface,  and  selected  structural 
features,  computed  using  the  TB  model  and  compared  with  DFT/LDA  results. 


TB 

DFT/LDA  a 

(100)  surface 

Ideal  surface 

r  (eV) 

2.13 

2.5 

1  X  1  relaxed  surface 

A  y  (eV) 

-0.02 

-0.03 

2  X  1  buckled  dimer  recons.  surf. 

A  y  (eV) 

-0.90 

-0.93 

dimer  bond  length 

d(  A) 

2.36 

2.19 

dimer  tilt  angle 

e 

17° 

15° 

(111)  surface 

Ideal  surface 

y  (eV) 

1.31 

1.39 

1  X  1  relaxed  surface 

A  y  (eV) 

-0.02 

2X2  T4  adatom  reconstruction 

A  y  (eV) 

-0.45 

-0.27  -  -0.30 

2X2  B2  adatom  reconstruction 

A  y  (eV) 

-0.04 

-0.10 

2X2  H3  adatom  reconstruction 

A  y  (eV) 

-0.12 

-0.16  -  -0.25 

‘From  Refs.  33-36  and  39. 
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FIG.  7.  Mean-squared  displacement  of  atoms  in  the  diamond 
lattice  as  a  function  of  temperature.  The  points  are  computed  from 
MD  simulations,  the  solid  line  is  a  linear  fit  going  through  the 
origin,  and  the  dashed  line  is  a  line  with  a  slope  corresponding  to 
the  experimental  measurement  as  discussed  in  the  text. 

bonds  with  its  nearest  neighbors.  No  DFT/LDA  results  for 
adatom-surface  bond  lengths  are  available.  The  rest  atom 
makes  three  2.40  A  long  bonds  with  angles  of  101.1°,  as 
compared  with  2.34  A  long  bonds  with  angles  of  104.9°  in 
DFT/LDA  calculations.  In  the  B2  site  the  TB  model  adatom 
makes  two  bonds  of  length  2.33  A  with  the  surface.  The 
binding  energy  in  this  configuration  is  underestimated  as 
compared  to  the  DFT/LDA  result.39 

D.  Finite  temperature  properties 

To  examine  some  finite  temperature  properties  of  the  dia¬ 
mond  structure  crystal  we  used  the  NRL  TB-MD  molecular- 
dynamics  package40  developed  by  Kirchhoff  to  evolve  a  512- 
atom  unit  cell  at  constant  energy  for  2000  molecular- 
dynamics  time  steps  (each  step  corresponds  to  2.0  fs), 
varying  the  initial  kinetic  energy  of  the  atoms.  From  the 
resulting  positions  we  computed  the  mean-squared  displace¬ 
ment.  This  quantity,  plotted  against  the  temperature  mea¬ 
sured  in  the  sample,  is  shown  in  Fig.  7.  A  linear  fit  of  the 
mean-squared  displacement  gives  a  slope  of  1.72 
X  105  A  2/K.  In  comparison,  the  Debye  temperature,  ex¬ 
tracted  from  experimental  measurements  of  the  temperature 
dependence  of  x-ray  diffraction  peak  broadening,41  corre¬ 
sponds  to  a  mean-squared  atomic  displacement  temperature 
coefficient  of  1.75X  10s  A  2/K.42 

From  the  Fourier  transforms  of  velocity  autocorrelation 
functions  calculated  during  MD  simulation  we  obtained  the 
vibrational  density-of-states.43  The  phonon-dispersion  curves 
were  extracted  from  the  Fourier  transform  of  the  velocity  and 
position  dependent  autocorrelation  function  for  a  given  po¬ 
larization  and  k  vector.43  The  vibrational  densities-of-states 
from  two  MD  simulations,  one  at  300  K  and  one  at  1500  K, 
are  plotted  in  Fig.  8.  The  overall  shapes  are  quite  similar, 
although  the  high-temperature  curve  is  smoother.  The  peaks 
are  shifted  to  lower  frequencies  at  1500  K,  indicating  a  soft¬ 
ening  of  the  vibrational  modes.  The  corresponding  phonon- 
dispersion  relations  are  plotted  in  Fig.  8,  and  compared 
against  experimental  values  extracted  from  the  graphs  of 
Ref.  44.  The  dispersion  relations  from  the  300  K  MD  simu¬ 
lation  are  in  good  agreement  with  experiment.  The  only  sig¬ 
nificant  errors  are  an  underestimate  of  the  frequencies  of  the 
second  branch  between  the  T  and  L  points  [A  , ( A ) ]  and  an 
overestimate  of  the  frequencies  of  the  highest  two  branches 
near  the  X  point  [A5(0),  A2'(0),  X2(0),  and  Xi(O)]. 


v(THz) 


FIG.  8.  Vibrational  densities-of-states  (upper  panel)  and  phonon 
dispersion  curves  (lower  panel)  extracted  from  the  velocity-velocity 
correlation  function  computed  during  a  molecular-dynamics  simu¬ 
lation  at  300  K  (solid  line)  and  1500  K  (dashed  line).  Dots  are 
experimental  data  measured  at  or  below  300  K,  extracted  from  the 
figures  in  Ref.  44. 

Some,  although  not  all,  of  the  vibrational  frequencies  in  the 
1500  K  MD  simulation  are  lower  than  at  300  K,  a  result  that 
is  consistent  with  the  differences  between  the  vibrational 
densities-of-states  in  Fig.  8.  The  high-frequency  branches 
have  the  largest  shifts,  with  many  of  them  shifting  down  by 
0.5  THz,  about  5%.  The  low-frequency  branches,  with  the 
exception  of  the  A  direction,  also  shift  down,  with  a  particu¬ 
larly  prominent  shift  of  the  second  branch  at  the  K  point 

[X3(A)]. 

E.  Amorphous  structures 

As  an  illustration  of  the  expanded  modeling  capabilities 
offered  by  the  present  TB  parametrizations  for  Si,  we  study 
the  electronic  properties  of  bulk  and  surface  structures  in 
amorphous  Si  (a-Si).  We  emphasize  at  the  outset  that  the 
following  discussion  is  concerned  mostly  with  demonstrating 
the  capabilities  of  the  approach,  rather  than  the  physics  of 
a-Si,  the  latter  being  a  broader  problem  beyond  the  scope  of 
the  present  paper. 

Based  on  many  experimental  and  theoretical  studies,45,46 
it  is  widely  accepted  that  a-Si  has  the  basic  structure  of  a 
continuous  random  network  (CRN)  of  tetrahedrally  bonded 
atoms,47-49  but  the  question  of  defects  has  been  the  subject 
of  considerable  debate  in  recent  years50,51  and  remains 
controversial.52  Experimental  results  appear  to  favor  under¬ 
coordinated,  three-fold  bonded  atoms  as  the  dominant  de¬ 
fects  (so-called  “dangling  bonds”).  On  the  other  hand,  the¬ 
oretical  simulations,  using  a  wide  variety  of  ab  initio , 52,53 
semiempirical54  56  and  empirical57-63  methods,  consistently 
produce  both  undercoordinated  as  well  as  overcoordinated 
(fivefold  bonded)  defects,  with  a  significant  preference  for 
the  latter.  The  type  of  bonding  arrangements  at  the  surface  of 
a-Si  is  even  less  clear  than  in  the  case  of  bulk  defects,  since 
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surface-specific  measurements  are  not  readily  available.  It 
has  been  reported  that  ab  initio  relaxation  of  a  bulk- 
terminated  CRN  model  produces  a  surface  with  roughly 
equal  numbers  of  threefold  and  fivefold  bonded  atoms.64  De¬ 
viations  from  the  tetrahedral  bonding  pattern,  either  in  the 
bulk  or  at  the  surface,  are  crucial  in  determining  the  elec¬ 
tronic  properties  of  the  material  because  they  introduce  states 
in  the  gap.52-54'64’65 

In  order  to  study  the  electronic  properties  of  bulk  and 
surface  a-Si,  we  first  prepare  a  32.9  A  X65.8  A 
X  16.5  A  bulk  amorphous  sample  with  1728  atoms.  Be¬ 
cause  such  a  large  sample  is  impractical  to  generate  using 
any  electronic  structure  based  simulation  method,  we  simu¬ 
late  the  quenching  of  the  liquid  with  an  interatomic 
potential,66  following  a  procedure  similar  to  the  one  used  in 
Ref.  62.  The  resulting  structure  has  over  96%  tetrahedral 
coordination,  with  only  fivefold  coordinated  defects. 

To  model  the  surfaces  of  a-Si,  we  considered  two  quali¬ 
tatively  different  1728  atom  slabs.  The  “cleaved  sample”  is 
created  from  the  bulk  structure  by  turning  off  the  periodic 
boundary  conditions  in  the  third  direction.  The  “quenched 
sample”  is  created  directly  from  the  liquid  phase  by  turning 
off  periodic  boundary  conditions  in  the  third  direction  and 
quenching  the  resulting  liquid  slab  with  the  interatomic  po¬ 
tential.  Not  surprisingly,  the  quenched  surfaces  are  slightly 
rougher  (by  about  1  A)  than  the  cleaved  surfaces.  The 
cleaved  surfaces  regions  contain  mostly  (65%)  fourfold  co¬ 
ordinated  atoms,  with  a  predominance  of  threefold  (29%) 
over  fivefold  (6%)  coordinated  atoms.  On  the  other  hand,  the 
quenched  surfaces  have  somewhat  higher  fourfold  coordina¬ 
tion  (72%),  with  many  (27%)  fivefold  and  almost  no  (1%) 
threefold  coordinated  atoms.  Note  that  these  surfaces  are 
fully  reconstructed,  and  therefore  ‘  ‘bulk’  ’  concepts  of  defects 
do  not  apply;  many  of  the  fourfold  coordinated  atoms  do  not 
have  tetrahedral  bond  angles,  and  the  fivefold  coordinated 
atoms  tend  to  appear  in  clusters  near  the  top  layer  of  the 
surface,  rather  than  as  isolated  floating  bonds  below  the  top 
layer. 

While  the  size  of  these  samples  makes  calculating  their 
electronic  structure  with  first-principles  methods  impractical, 
the  sp^d5  basis  TB  parametrization  described  in  Sec.  II 
makes  such  a  study  feasible.  In  Fig.  9  we  compare  the  elec¬ 
tronic  DOS  of  the  three  amorphous  samples  computed  with 
the  TB  model  (see  also  Fig.  3  for  the  diamond  structure 
crystal  electronic  DOS).  The  bulk  amorphous  sample  DOS  is 
much  smoother  than  the  crystal,  in  agreement  with  other 
simulation  results.53'67  There  are  only  two  peaks  in  the  va¬ 
lence  band,  corresponding  to  the  highest  and  lowest  peaks  of 
the  DOS  of  the  crystal.  Despite  the  structural  differences 
between  the  two  surface  models,  their  overall  DOS  is  quite 
similar,  showing  that  the  electronic  signatures  of  undercoor¬ 
dinated  and  overcoordinated  atoms  at  the  amorphous  surface 
are  difficult  to  distinguish.  This  result  is  consistent  with  the 
arguments  of  Pantelides  for  bulk  defects.50  The  DOS  of  the 
surface  samples  differ  from  that  of  the  bulk  sample  essen¬ 
tially  only  in  the  gap  region.  This  indicates  that  all  surface- 
related  defects  produce  gap  states,  consistent  with  analogous 
results  for  bulk  defects. 52,54,65  It  is  also  interesting  that  the 
surface  DOS  above  the  gap  region  is  depleted  relative  to  the 
bulk  one.  A  more  detailed  analysis  of  these  results  will  be 
given  elsewhere.68  Here,  we  wish  to  point  out  only  the  effi- 
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FIG.  9.  The  electronic  TB  DOS  computed  with  the  sp3d5  pa¬ 
rametrization  for  the  bulk  and  surface  a-Si  samples.  The  lower 
curve  (dashed  line)  is  the  DOS  for  the  bulk  amorphous  sample.  The 
two  upper  curves  show  the  excess  DOS  associated  with  the  sur¬ 
faces,  plotted  as  the  difference  between  the  cleaved  sample  DOS 
and  bulk  amorphous  DOS  (dotted  line),  and  the  difference  between 
the  quenched  sample  DOS  and  the  bulk  amorphous  DOS  (solid 
line). 

ciency  of  the  approach  in  generating  reliable  electronic  struc¬ 
ture  information  for  large  systems. 

IV.  SUMMARY 

We  have  applied  the  NRL-TB  method  to  generate  TB 
models  for  Si  that  were  fit  to  LAPW  results  of  a  small  num¬ 
ber  of  high-symmetry  crystal  structures.  We  found  that  the 
resulting  Hamiltonians  are  transferable  to  a  much  wider 
range  of  geometries.  A  model  with  a  nonorthogonal  yp3  ba¬ 
sis  reproduces  DFT/LDA  and  experimental  measurements 
for  a  wide  range  of  material  properties,  including  elastic  con¬ 
stants  and  phonon  frequencies,  point  defect  formation  ener¬ 
gies,  and  surface  energies  and  reconstructions.  In  fact,  this 
TB  model  is  as  good  or  better  at  describing  the  energetics  of 
point  defects  than  some  models  that  included  such  structures 
in  their  fitting  process.  It  is  also  the  only  nonorthogonal  TB 
model  we  are  aware  of  that  correctly  describes  the  energy 
sequence  of  different  adatom  configurations  on  the  (111)  sur¬ 
face  of  silicon.  The  ability  of  the  model  to  accurately  de¬ 
scribe  such  diverse  systems,  despite  having  been  fitted  to  a 
small  number  of  high-symmetry  crystal  structures,  increases 
our  confidence  that  the  model  captures  the  essential  physics 
of  bonding  in  solid-state  silicon  systems. 

The  efficiency  of  the  model  makes  it  possible  to  study 
finite  temperature  properties  of  large  silicon  systems  through 
molecular-dynamics  simulation,  an  application  that  would  be 
impractical  with  DFT/LDA  methods.  We  have  shown  that 
the  model  reproduces  experimental  results  for  atomic  mean- 
squared  displacements  as  a  function  of  temperature  in  bulk 
silicon.  Phonon  densities-of-states  and  dispersion  curves  ex¬ 
tracted  from  MD  simulations  at  different  temperatures  show 
good  agreement  with  experiment  at  low  temperatures,  and  a 
substantial  softening  of  many  modes  at  higher  temperatures. 
By  adding  d  orbitals  and  modifying  the  fitting  data  set,  we 
obtained  a  model  that  accurately  reproduces  both  the  valence 
and  conduction  bands  of  silicon  in  the  diamond  structure,  at 
the  price  of  deterioration  in  the  accuracy  of  the  energetics. 
This  sp3d5  parametrization  makes  possible  the  study  of  the 
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electronic  structure  of  amorphous  systems  with  nearly  2000 
atoms. 
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